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INTRODUCTION 


Many  problems  in  science  and  engineering  require  the  solution  of  a 
system  of  partial  differential  equations  on  a  complex  two-,  or 
three-dimensional  domain.  In  order  to  solve  such  a  problem 
numerically,  a  boundary -fitted  coordinate  system  is  usually 
constructed.  For  geometries  of  practical  importance,  it'  is  often 
necessary  to  partition  the  domain  and  to  fit  a  separate  coordinate 
system  to  each  subregion.  .Patching  these  subsystems  together  to  form  a 
composite  system  can  produce  irregular  grid  structures  where  the 
numerical  scheme  of  solution  must  be  modified. 

7 

Methods  for  treating  nonstandard  (‘"special")  grid  points  in 
two-dimensional  conf igurations  were  developed  and  reported  earlier  in 

the  course  of  this  project.  -{See  Ref. . 1~3)'.  Verification  of  these 

methods  through  computer  experimentation  and  extension  of  the  methods' 

to  special  points  in  three-dimensional  grids  were  the  primary 

/'  -  ‘ 

objectives  of  the  project  reported  here.  — ■  '  '  J  •  1 


PROJECT  RESULTS 


The  same  general  approach  us.ed  to  develop  and  analyze  methods  for 
treatment  of  special  two-dimensional  grid  structures  was  applied 'to 
three-dimensional  configurations.  The  characterization  of  special 
points  in  three'  dimensions  is  a  direct  and  straightforward  extension  of 
that  reported  for  two-dimensional  geometries  (Ref.  1).  Ordinary 
three-dimensional  grid  structures  are  depicted  in  Fig.  1  and  2.  Such 
structures  may  not  exist  at  some  boundary  points  or  where  two  or  more 
subregions  of  a  composite  grid  meet.  Examples  are  shown  in  Fig.  3  and 
4.  These  nonstandard  grid  structures  require  an  adjustment  of 
numerical  solution  schemes  whicn  are  based  on  the  usual  grid-point 
orientation  shown  in  Fig.  1  and  2. 

As  in  the  case  of  two-dimensional  geometries,  a  finite-difference 
scheme  can  be  applied  at  a  special  point  in  three  dimensions  if  the 
scheme  is  rephrased  in  terms  of  a  suitable  local  coordinate  system. 
This  local  system  must  have  an  ordinary  structure  (as  in  Fig.  1  or.  2) 
and  should  not  be  too  stretched  or  skewed.  The  appropriate  local 
coordinate  systems  for  the  special  points  in  Fig.  3  and  4  .are 
illustrated  in  Fig.  5  and  6.  It  should  be  noted  that  for  some 
configurations,  a  satisfactory  local  coordinate  system  may  not  exist. 
An  .example  is  shown  in  Fig.  7.  In  such  cases,  a  f ini te- volume  approach 
must  be  used  to  approximate  derivatives'  at  the  special  point. 

Irregular  grid  structures  can  generally  be  treated  in  a  simple  and 
straightforward  manner  when  a  finite-volume  approach  is  used.  Then,  is 
is  the  configuration  of  grid  cells  which  is  at  issue.  Extension  of  the 
two-dimensional  analysis  of  Ref.  1-3  shows  a  special  cell  in  three 


dimensions  to  have  other  than  six  faces  or  a  vertex  common  to  'other 
than  eight  cells.  The  latter  type  of  special  cell  requires  special 
attention  only  when  a  derivative  must  be  represented,  on  a  cell  face. 

In  direct  analogy  to  the  approach  for  two  dimensions  detailed  in 
Ref.  2  and  3  (attached  herewith),  special  cells  are  treated  by  imposing 
the  integral  form  of  the  governing  equations  in  physical  coordinates, 
rather  than  in  terms  of  the  transformed  (curvilinear')  coordinates, 
within  the  special  volume  element.  Cell-centered  values  of  the  unknown 
functions  are  sought.  A  derivative  within  the  cell  is  approximated  by 
a  -surface  integral  over  the  cell  boundary  which  is  equivalent  to  the 
average  of  that  derivative  over  the  cell.  In  particular,  identity  (1) 
is  used: 


7 f  -  y-  I  Vf  dc  «  yr  f  ■  f  n  da, 

~  ay3  '  -!D  ~  V  ‘  3D  ~ 


where  7  is  the  volume  of  cell  D,  dc  is. an  increment  of  volume,'  5D  is 
the  surface  bounding  ceil  D,  and  n  is  the  unit  outward  normal  to 
surface-area  element  da.  The  value  of  the  unknown  function,  f,  on  a 
cell'  face  is  equated  to  an  average  of  the  function's  values  at  the 
centers  of  the  cells'  sharing  that  face.  The  value  of  a  derivative, 
e.g.  — — ,  cn  a  cell  face  is  equated  to  the  average  of  tr.e  derivative's 
values '"at  the  vertices  of  that  ceil  face.  These  vertex  values  are 
obtained  by  appivir.g  identity  (1)  to  a.n  auxiliary  cell  centered  at  each 
vertex.  The  cumbersome  formulas  for  these  expressions  in  three 


dimensions  are  c 


mitted  here,  but  they  can  be  derived  by  directly 


extending  the  formulas  for  'two-dimensional  cells  given  in  Ref.  2. 
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Abstract 

Complex  physical  reqions  often  necessitate 
the  use  of  several  rectangular  computational 
regions.  The  corresponding  physical  subregions 
may  be  joined  along  their  boundary  or  overlapped. 
Techniques  are  described  for  the  numerical 
solution  of  partial  differential  equations  on 
either  type  of  composite  grid  system. 

Introduction 

Problems  in  constructing  a  computational 
grid  often  arise  due  to  the  complex  geometric 
configurations  encountered  in  practical  fluid 
flow  simulations.  For  example,  an  aircraft 
is  composed  of  many  individual  components  such 
as  fuselage,  wings ,  tail,  fins,  pylons,  and 
nacelles.  Although  each  component  may  be 
geometrically  simple,  and  a  grid  about  the  in¬ 
dividual  components  may  be  easy  to  construct, 
constructing  a  grid  about  the  complete  con¬ 
figuration  is  more  difficult.  Another  example 
in  which  complex  regions  are  encountered  may 
be  found  by  examining  the  design  of  current 
turbomachinery.  The  . problem  of  constructing  a 
grid  is  compounded  by  the  accuracy  requirements 
of  the  numerical  algorithm.  In  order  to  control 
trie  truncation  error  of  the  algorithm,  grid 
points  must  be  concentrated  near  boundary  layers, 
shocks,  and  other  regions  where  solution  gradi¬ 
ents  are  large.  It  is  thus  apparent  that 
where  a  boundary- fitted  curvilinear  coordinate 
system  is  employed  in  the  solution  of  many 
problems  of  practical  interest,  several  rec¬ 
tangular  computational  regions  -mst  be  used. 

The  grid  which  is  defined  by  this  merging  of 
several  coordinate  systems  is  called  a  composite 
grid. 

Two'  approaches  have  developed  in  construct¬ 
ing  composite  grid  systems.  The  first,  and 
still  the  most  popular,  is  to  have  the  grid 
lines  pass  smoothly  from  one  subregion  into  the 
next.  ',2  The  well-kncwn  chain  rule  formulas 
can  be  used  to  transform  from  the  physical  to 
the  computational  variables  at  the  interior 
points  of  each  subregion  as  well  as  at  points 


J 

lying  along  the  sides  of  two  subregions.  However 
the  vertices  present  a  problem  whenever  the 
Jacobian  of  the  transformation  is  zero  or  is 
undefined.  This  same  problem  may  also  occur  at 
boundary  points  of  the  physical  region.  One  of 
the  objectives  of  this  report  is  to  classify  the 
types  of  singular  or  unusual  points  which  may 
occur  and  indicate  how  to  derive  consistent 
difference  approximations  at  such  points.  The 
second  approach  in  working  with  composite  grids 
Is  to  overlap  the  grids  for  the  physical  sub- 
regions.  In  this  case  the  solution  values  are 
transferred  from  the  Interior  points  of  one  grid 
to  the  boundary  points  of  another  grid  by  an 
interpolation  formula.  Many  different  inter¬ 
polation  schemes  are  available  and  their  effect 
on  the  error  In  the  numerical  solution  may  be 
significant.  Overlapping  grids  have  been  gaining 
in  popularity  due  to  the  Inherent  freedom  in 
generating  each  subregion  grid  without  considera¬ 
tion  of  the  remaining  grids. 

Computational  techniques  and  comparative 
results  on  composite  grids  will  be  presented  in 
the  solution  of  simple  two-dimensional  model 
problems.  In  the  case  where  grid  lines  con¬ 
tinue  smoothly  between  subregions,  a  finite- 
volume  method  Is  compared  with  sn  ADI  scheme 
using  the  usual  finite-difference  approximations. 
The  first  method  illustrates  the  use  of  finite- 
volume  formulations  for  solving  partial  differen¬ 
tial  equations  with  second-order  derivatives, 
while  the  second  demonstrates  an  imp!er«ntation 
of  ADI  schemes  on  composite  grids  with  singular 
points.  Solutions  are  also  presented  from  com¬ 
putations  on  overlapping  grids.  The  questions 
to  be  discussed  there  are:  What  effect  does  the 
Interpolation,  procedure  have  on  the  accuracy  of 
the  numerical  solution?  How  does  the  accuracy 
vary  with  the  choice  of  the  particular  inter¬ 
polation  formula?  Specifically,  comparisons  are 
presented  using  bivariate  interpolation  on 
quadrilaterals  3,4t  Taylor  polynomials  5,  and  the 
triangular  interpolation  schemes  corr-only  used 
with  finite  element  methods. 


Special  Points  In  Continuous  Grids 


•This  work  was  supported  by  the  Army  Research 
Office  under  Grant  DAAG  29-83-K-0101  and  by 
NASA  Langley  Research  Center  under  Grant 
NSG-1577. 


Field  points  which  require  special  attention 
frequently  arise  when  numerically  solving  flew 
problems  on  a  composite  curvilinear  coordinate 
System  consisting  of  continuously  joined  segments. 
This  is  true  whether  the  numerical  method  of 
solution  involves  a  finite-volume  formulation  or 
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a  finite-difference  scheme.  The  special  points 
encountered  in  the  two  approaches  are,  however, 
different  in  nature,  hence  the  two  cases  are 
addressed  separately  in  the  following  discussion. 

In  a  finite-volume  scheme  on  a  two-dimen¬ 
sional  composite  grid,  the  unusual  points  which 
arise  are  centers  of  Interior  grid  cells 
having  other  than  four  edges  or  having  a  vertex 
shared  by  other  than  four  cells.  Examples  are 
illustrated  below. 


The  first  type  of  cell  can  result  when  sub- 
regions  of  the  segmented  domain  are  Joined 
between  grid  lines;  the  second  can  occur  when 
segments  are  joined  along  grid  lines,  but  this 
cell  requires  special  care  only  when  second 
derivatives  are  approximated  therein.  In 
either  case,  treatment  of  the  associated 
"special  points"  consists  simply  of  replacing 
the  usual  difference  representations  associated 
with  four-sided  cells  having  vertices  common 
to  three  other  cells  by  expressions  derived  for 
a  general  fi-sided  cell. 

Recall  that  a  flni te-volume  formulation 
solves  for  cell-centered  function  .values  and 
approximates  derivatives  at  a  cell  center  by 
line  Integrals  about  the  cell  boundary  which 
are  equivalent  to  averages  over  the  cell.  In 
particular,  identity  (1)  is  used. 


/'  5.  ds 


For  an  N-sided  cell  of  area  A  with  cartesian 
centroid  P  *  (p,.  P,).  vertices  V’  *  1 

(vi,  4)  1*  1  2,  . . .  •  N,  ard  edges  s' 

'  1  2  ,-*1  N+i  i 

joining  V  and  V  (V  «  V  )  along  which 
a  function  f  and  its  first  partial  derivatives 
are  constant,  this  approach  gives 


<  *  A'1  i  f*1  (v’*1-  vj) 
1*1  ‘  L 

s-5  - 

1*1 

4  .  A fS\v}  -  v}+1) 
yy  fTj  y  1  1 


where  the  superscripts  on  f  and  its  derivatives 
Indicate  the  point  or  face  of  evaluation. 

Simple  forward/backward  differencing  can  le 
us^d  to  approxjmate  edge  values 

fs,  fs,  and  fs  when  working  on  a  grid  whose 
x  y 

Interior  Is  "rectangular"  In  nature  In  that 
right/left  and  up/down  orientations  can,  be 
followed  through  the  grid  In  an  unambiguous  way. 
Such  a  technique  was  used  in  the  numerical  ex¬ 
ample  discussed  later.  On  a  more  general  grid, 
an  obvious  way  to  approximate 

fs  is  to  average  the  center  values  of  the  two 
cells  sharing  edge  $1.  This  same  averaging  scheme 
cannot  be  repeated  to  approximate 

1  4 

fs  ,  and  fs  ,  however,  without  rejecting  the 
x  y 

customary  strategy  of  avoiding  use  of  values  at 
points  which  are  not  Immediate  neighbors  of  the 
point  at  which  a  quantity  is  being  evaluated. 
Instead,  we  propose  the  averaging  technique: 


*5  ■  t  *f; 


i  i  vi  vi*i 

fy-?(fy+fy 


where  the  vertex  values  are  obtained  by  applying 
Identity  (1)  to  auxiliary  cells  formed  by 
joining  the  midpoints  of  the  edges  of  each  cell 
to  the  ceil  center.  To  nuke  this  .more  precise, 
let  V  be  a  vertex-  common  to  Q  cells  and  label 
the  cell  faces  emanating  from  V  as  k'  with 
midpoints 

H1  -  (m{.  m*)  1*1,  2  Q-  Then,  if  P  • 

(pj,  p*)  Is  the  center  of  the  cell  having  edges 

k1  and  k1+1,  and  If  f  is’take"  to  equal  fP  jln"3 
mV  and  pV+1  ,  the  first  partial  derivatives 
of  f  at  V  may  be  approximated  by 


■  a-’£  /(.;*'  -  -;> 


fV  .  *'4  *  mj*’)  , 


where  A  Is  the  area  of  the  2Q- faced  auxiliary 

cell  mVmV  ...  mW  Indicated  In  tne 
following  diagram. 


(2) 


The  usual  transformation  relations  applied  to 
the  local  coordinates  (l,  n)  at  a  special  point 
P  provide  accurate  difference  approximations  of 
derivatives  with  respect  to  the  physical  variables 
there.  This  technique  was  applied  to  the  special 
boundary  point  shown  in  figure  2  when  the 
associated  partial  differential  equatlon'was 
numerically  solved  using  a  finite-difference 
scheme. 


Expressions  (3),  (4),  and  the  average 


1  p 

2 


1-1 


/) 


mate  the  difference  formulas  In  12}  reasonable 
approximations  of  the  first  and  second  partial 
derivatives  of  a  function  f  at  the  center  of 
an  arbitrary  grid  cell.  Because  this 
technique  Is  generally  applicable,  the 
finite-volume  approach  Is  deemed  the  roost 
straight-forward  when  special  points  are  pres¬ 
ent  in  the  field.  In  the  alternative  finite- 
difference  approach.  It  Is  found  that  special 
points  which  arise  must  be  Individually  ex¬ 
amined  before  a  suitable  treatment  of  each 
can  be  identified.  Thus  no  single  set  of 
finite-difference  formulas  can  always  be 
applied  to  every  such  point  independent  of  its 
context,  as  will  now  be  discussed. 


The  special  points  which  are  encountered 
in  a  finite-difference  formulation  may  be 
recognized  as  those  grid  points  (cell  vertices) 
having  a  nonstandard  number  of  Imediate 
neighbors.  On  the  boundary  of  a  twc-cijnensional 
domain,  such  points  have  other  than  five 
neighbors  but  require  special  treatment  only 
if  Neumann  •'uundary  conditions  are  Imposed. 
Interior  grid  points  are  "special"  If  they  have 
other  th.3n  eight  irrediate  neighbors.  These 
can  occur  on  the  interface  between  segments 
of  a  composite  grid  which  are  joined  along 
grid  lines,  or  they  can  be  vertices  of  a  cell 
intersected  by  an  interface  in  the  caso  of 
subregions  joined  between  grid  lines  (see 
below) . 


In  all  cases,  the  usual  finite-difference 
approach  can  be  followed  if  the  difference 
for.mjlas  are  rephrased  in  terms  of  suitable 
local  coordinates  which  Identify  and  label  , 
the  standard  number  of  neighbors.  Thompson® 
has  described  such  local  systems  for  special 
boundary  points;  choices  appropriate  to 
various  interior  special  points  are  indicated 
in  Figure  1. 


Overlapping  Grids, 

When  overlapping  grids  are  constructed, 
any  boundary  point  of  a  component  grid  G,  which 
Is  not  a  physical  boundary  point,  must  lie  in 
a  cell  of  some  grid  H.  Regardless  of  the  differ¬ 
ence  equations  or  the  algorithm  used  to  solve 
the  difference  equations,  there  roust  be  a  trans¬ 
mission  of  Information  between  the  v.arious  grid 
systems.  Suppose  that  solution  values  on  H  are 
transmitted  to  the  grid  G  by  a  general  inter¬ 
polation  formula 

f(rQ)  •  £  a/tSj)  ,  (S) 

J*1 

where  r.  is  a  boundary  point  of  G  and  the  s.  are 
points  of  H.  At  least  some  of  the  s.  must  J 
be  Interior  points  of  1  and,  de-ending  on  the 
algorithm,  one  may  wish  to  require  a  sufficiently 
large  overlapping  so  that  all  the  s.  are  interior 
points.  For  example,  the  latter  req ui rement 
would  eliminate  possible  stability  problems  when 
solving  parabolic  equations.  When  the  value 
of  f  at  r0  is  replaced  in  the  difference 
equation  at  a  neighboring  point  by  th.e  inter¬ 
polated  value,  a  new  difference  equation  results 
which  will  have  a  different  local  truncation 
error.  It  can  be  shown  that  a  sufficient 
condition  for  consistency  in  tit  approximation 
of  the  first-order  derivatives  with  njspect 
to  the  physical  variables  x  and  y  is  that 


k  k  k 

£  .J  -  l.E  Yi  *  V  £  *jyj  *  yo  •  (6) 

J-l  J-l  j»0. 


where  Xg  and  y^  denote  the  coordinates.of  r^,  and 

x.  and  y.  denote  the  coordinates  of  s..  The 
J  3  — . 

difference  approximations  will  be  second-order 
accurate  if  in  addition 


(7) 


£ 

J*1 


Vj 


a.*  .y . 
j  yj 


v,,£ 


x.v.,^  a.y^ 
^ u  j-l  J  J 
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There  Is  a  loss  of  one  order  In  the  formal 
accuracy  of  the  second-order  derivative 
approximations. 

Several  different  interpolation  schemes  will 
now  be  reviewed.  The  first  scheme  Is  based  on 
the  approximation  of  a  Taylor  polynomial.  For 
each  boundary  point  rn  of  G,  find  the  closest 
point  s,  of  H.  The  values  of  the  function 

f  at  th4  four  neighbors  of  s. ,  say  s. .  s,, 

can  be  used  to  calculate  approximations  5 
of  the  first  derivatives  of  the  actual  solution. 
These  approximations  can  thus  be  substituted 
for  the  derivatives  In  a  linear  Taylor  polynomial 
at  s^ .  If  the  resulting  interpolation  formula 
is  written  in  the  form  of  (5),  then  the  coef¬ 
ficients  will  satisfy  (6).  If  the  original 
algorithm  is  at  least  second  order  accurate, 
then  there  is  no  loss  of  accuracy  In  replacing 
the  derivatives  by  differences.  Therefore, 
the  use  of  this  interpolation  formula  would 
give  consistent  approximations  of  first 
derivatives.  Clearly,  one  could  also  construct 
an  approximation  of  a  quadratic  Taylor  polynomial. 
However,  the  consistency  of  the  second  deriva¬ 
tive  approximations  cannot  be  guaranteed  unless 
the  local  truncation  error  in  the  tinite- 
difference  equations  exceeds  two.  Therefore,  if 
a  second-order  algorithm  is  being  implemented, 
it  is  not  clear  that  a  quadratic  polynomial  would 
give  better  results  than  a  linear  polynomial. 

The  next  interpolation  formula  does  not 
require  the  qu.-nerical  approximation  of  deriva¬ 
tives.  Let  rQ  belong  to  a  grid  cell  C  of  H 
with  vertices  us.,...,  s..  There  exists  a 
unique  bilinear  mapping^of  the  unit  square  onto 
the  cell  C.  The  point  rQ  must  be  the  image  of 
sene  point  in  the  unit  square  so  that 


rQI(l-u)'l-v)51  +  u(l-v)s2  ♦  v(l-u}s4  ♦  uvs3 

for  0  <_  u,  v  <_  1.  Since  the  coordinates  of  all 
points  are  xnown,  this  system  can  be  solved 
explicitly  for  u  and  v.  If  f  Is  also  assumed 
to  be  a  bilinear  function  of  u  and  v,  then 

• 

f(r0/*O-u)(l-v)f(Sj)  ♦■u(l-v)f(s2)  +  v(l-u)f(s.; 
+  uvf(s3).  (8) 

This  interpolation  formula  satisfies  (6).  This 
scheme  can  also  be  generalised  to  construct 
interpolation  formulas  of  arbitrary  order.  The 
procedure  for  biquadratic  interpolation  will 
be  briefly  described.  Let  Q  be  the  union  of 
four  grid  cells  containing  rQ.  There  exists  a 
unique  biquadratic  mapping  of  the  unit  square 
Such  that  the  points  (a,  v)  where  u*0,  1/2,  ■ 

1,  v=0,  1/2,  1,  map  to  the  vertices  of  Q.  Now 
the  coordinates  of  the  point  in  the  unit  square 
which  has  r  as  its  image  can  be  found  by  solving 
a  system  of°two  biquadratic  equations. 
Unfortunately,  these  equations  would  be  difficult 
to  solve  directly.  In  the  examples  which 
follow,  these  equations  were  solved  numerically 
using  Ne-tcn's  iteration  method.  Once  the 


approximate  values  of  u  and  v  have  hee  calcu¬ 
lated,  a  biquadratic  interpolanc  of  the  following 
form  is  computed. 

9 

f(r0 )  *  £  Qj  (u.  v)f(Sj )  (9). 

J*1  , 

Here  q.  Is  the  biquadratic  Lagrange  interpolating 
polynomial  which  assumes  the  value  1  at  the 
point  of  the  unit  sq-are  which  maps  to  s,  and' 
vanishes  at  the  points  which  have  other  grid 
points  of  Q  as  images.  This  interpolation  for 
mula  satisfies  both  (6)  and  (7).  Thus,  improved 
results  would  be  anticipated  using  biquadratic 
interpolation.  It  should  be  noted  that  the 
above  interpolating  functions  are  bilinear  and 
biquadratic  functions  of  the  variables  u  and  v 
which  are  sometimes  referred  to  as  Isoparametric 
coordinates.  Bilinear  and  higher-order  Inter¬ 
polation  on  the  physical  variables  x  and  y  is 
generally  not  feasible  unless  the  cell  sides  are 
aligned  with  the  coordinate  axes. 

Partly  due  to  the  need  for  an  iterative 
method  in  determining  the  interpolation 
coefficients,. it  was  decided  that  interpolation 
on  triangular  regions  would  be  investigated. 

Such  interpolation  is  popular  In  finite-element 
analysis,  and  it  can  be  readily  adapted  to  the 
present  problem.  If  rn  belongs  to  a  grid  cell 
C,  then  either  diagonal  of  (1  determines  a  tri¬ 
angular  region  containing  tq.  Select  one  of  these 
triangular  regions  and  let  its  vertices  be  de¬ 
noted  as  s. ,  s?,  S-.  In  this  case,  the  three 
equations  in  ^(6) determine  the  coefficients 
a.  of  the  linear  i..-er?olating  formula  (5). 

JThis  is  referred  to  as  linear  Interpol ation 
since  it  Is  equivalent  to  Interpolating  by  a 
linear  function  of  x  and  y  with  values  given 
at  s, ,  s?,  and  s,.  A  quadratic  interpolation 
can  5e  constructed  from  three  points  at  the 
vertices  of  a  triangular  region  and  three  points 
along  the  sides.1  Such  a  configuration  arises 
if  the  grid  points  of  H  are  selected  as  in  this 
figure. 


Again  the  coefficients  of  the  interpolation 
formula  are  solutions  of  a  linear  system  of 
equations,  namely,  the  equations  in  (6}  and  (7). 
Then  the  i -.terpol  a  ted  value  is  obtained  by  using 
these  coefficients  in  (5). 

So  far  only  accuracy  considerations  have 
been  investigated.  When  solving  elliptic 
and  parabolic  equations.  Iterative  convergence 
and  stability  are  also  of  interest.  If  only  the 
Interpolation  formula  (5)  is  considered,  con¬ 
vergence  and  stability  would  be  guaranteed  if 
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the  cueff icivnts  are  non-negative.  Tnis  wools 
therefore  be  a  desirable  condition  although 
certainly  not  necessary.  The  linear  arc  bi- 
’irear  are  the  only  interpolation.  formulas 
with  nor-negati  .e  coefficients..  It  car.  also 
be  noted  that  the  convergence  or  stability  of 
the  finite-difference  algorithm,  would  be  least 
affected  if  there  were  sufficient  overlap  of 
grids  so  that  all  the  s  are  interior  grid  points 
of  H.  This  is  because  all  the  f ( s . )  values 
would  be  updated  when  the  Interpolation  procedure 
was  implemented; 


Con-y u tatlor.ai  ttamples 

The  prcble~s  considered  here  are,  for  the 
most  port,  simple  model  problems  where  an 
enact  solution  is  known.  Hence,  the  accuracy 
of  the  various  techniques  discussed  in  .the  pre¬ 
vious  sections  can  be  compered.  The  physical 
regions  are'also  quite  simple,  but  they  do 
exhibit  the  computational  problems  Inherent  In 
more  complex  conf 1 gura t ions . 

The  advantages  of  the  cel  1 -centered  finite 
velum*  method  in  the  solution  of  problems  on 
regions  with  special  points  has  already  been 
discussed,  recent  literature  cortalns  many 
applications  of  the  method  to  the  solution  of 
hgp-rtclic  conservation  laws.  In  the1  case  of 
fo-r-sided  ceils,  the  method  can  a’se  be 
applied  In  a  straightforward  manner  to  the 
solution  of  parabolic  ana  elliptic  equations 
by  using  the  for-ard/back-ard  differences  which 
were  mentioned  earlier.  A  second-on den  ap¬ 
proximation  is  .obtained  by  Switching  tr.e 
direction,  of  differencing  In  the  't-o  difference 
cceratip-s  which  must  be  applied  to  ce-e-ate  an 
approximation  of  a  second- -order  derivative. 

In  the  results  of  this  report,  the  four 
possible  difference  quotients  obtained  by  all 
com- 'rations  of  forward  and  baci-srd  differencing 
have  been  averaged  to  yield  a  second-order 
d; ffer ence  expression  which  would  also  reflects 
any  symmetry  in  the  grid.  It  was  observed 
that  considerable  simplification  in  imposing 
boundary  conditions  resulted  if  only  two  of  the 
difference  expressions  were  averaged  at  the' 
bp.rpary.  If  Se-nvar-n  boundary  conditions  -ere 
specified,  differencing  in  the  direction  ex¬ 
terior  to  the  region  was  performed  only  in  the 
second  appl ication  of  the  difference  operator  and 
net  1r,  the  computation  of  .the  first  derivative. 

U"  en  Cirich’et  boundary  conditions  were  Inpcseo, 
differencing  in  the  exterior  direction  was  per- 
fer-xed  only  on  the  first  application  of  the 
di f ference  operator.  This  converlence  eliminates 
the  need  for  false  boundary  points  ar-d  permits 
direc'  substitution  of  either  Dirich', et'  or 
Neumann  data  into  the  difference  equations. 

The  finite  volume  method  has  been  used 
to  solve  Laplace's  equation  in  a  region  about 
an  ellipse.  The  grid  which  was  used  is  shown 
ir  Figure  2.  A  singularity  occurs  at  a  vertex 
of  the  ellipse  where  three  gr’d  lines  intersect 
at  a  single  point.  In  order  to  com, pare  with  a 
known  solution,  the  boundary  values  of  the 
potential  function  for  ideal  flc-w  about  an 


elliptic  cylinder  were  i -posed.'  Specifically, 
the  boundary  value  prosier 


p  +  p  •  0  in  the  interior 
rxx  ryy 

Pn  ■  0  on  the  ellipse  (1 0} 

p  ■  x  at  the  outer  boundary 


was  solved  for  the  region  ir.  Fig-re  2.  Note 
that  only  half-cells  are  depicted  at  the  outer 
boundary  where  Cirichlet  'boundary  conditions  are 
imposed  at  cell  centers.  The  variation  of 
the  computed  surface  potential  along  the  ellipse 
was  In  agreement  with  the  theoretical  results. 
However,  the  solution  appeared  to  be  translated 
by  a  small  constant  term.  This  was  traced 
to  a  larger  thanexpected  error  In  the  coarse 
grid  region  near  the  outer  boundary.  The  error 
at  the  outer  boundary  was  subtracted  from 
the  numerical  solution  and  the  resulting  surface 
potential  was  compared  with  the  theoretical 
solution.  The  absolute  error  In  the  surface 
potential  Is  plotted  in  Figure  3  beginning  with 
tne  vertex  at  x*C  and  ending  at  the  vertex 
x=l  where  the  singularity  occurs.  For  the 
number  of  grid  points  which  were  used,  one  might 
expect  more  accurate  results.  This  error  was 
partially  due  to  the  large  error  at  the  Outer 
boundary  which  needs  further  investigation. 

The  computed  solution  was  also  not  defined  on' 
the  surface  and  thus  a  linear  extrapolation 
of  .interior  values  was  used  to  obtain  the 
numerical  values  which  were,  compared  with  the 
act-al  solution.  A  further  source  of  error 
was  the  fact  that  the.  free  stream,  boundary 
condition  was  applied  only  one  and  a  half  chord 
lengths  from  the  ellipse. 

The  computation  of  the  potential  function 
for  flew  about  the  same  elliptic  cylinder  was 
also  carried. out  using  the  traditional  finite- 
difference  method.  Point  SOP  iterative  methods 
are  frequently  used  for  solving  elliptic 
equations  on  composite  grids.  That  method  was 
used  in  the  above  finite -volume  calculations 
and  could  also  be  used  here.  However,  It  was 
decided  Instead  to  use  an  ADI  method.  This 
example  Illustrates  the  special  treatment  needed 
at  singular  points  3nd  the  successful  application 
of  implicit  updating  along  boundary  cjrves  in 
rectangular  subregions.  In  this  case  there  is 
only  one  computational  rectangle,  but  there 
are  two  re-entrant  boundary  segments  which  map 
to  the  same  interior  grid  line. 

Laplace's  equation,  as  well  as  other 
elliptic  equations,  can  be  transformed  to  com¬ 
putational  variables  and  written  in  the  form 

SPU  +  b  Ptn  +  C  Pnn  +  d  Pt  +  6  Pn  '  °  (11) 
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The  -31  method  which  was  used  to  solve  this 
eolation  is  the  Douglas-Fachfprd  algorithm 
-r'ch,  after  derivatives  are  replaced  by 


fe-ences,  can  be  written  as 
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C*  course  these  equations  were  appropriately 
modified  at  the  boundary  to  reflect  the 
i.e.-ann  boundary  conditions  in  (10).  Suppose 
that 'the  computational  region  is  the  unit  square 
nth  the  image  of  the  ellipse  lying  along  the 
C’O  coordinate  line.  The  grid  is  similar 
to  Figure  2.  Since  Dirichlet  boundary  conditions 
are  imposed  at  grid  points,  no  half  cells  are 
used  at  the  outer  boundary.  The  computations 
in  (IT)  and  (13)  can  be  carried  out  in  the 
following  steps: 

(i)  Calculate  p*  on  the  ellipse. 

(ii)  Using  the  value  of  p*  at  the  singular 
point,  calculate  p*  on  the  re-entrant 
segment  of  the  coordinate  line. 

(iii)  Calculate  p*  on  the  remaining  i= 
constant  coordinate  lines. 

(iv)  Calculate  p^1*'^  on  all  n ’constant 
coordinate  lines  whicn  intersect  the 
ellipse  at  any  point  other  than  the 
,  s ;  -  g'ul  ar  pci,nt. 

(v)  Calculate  p' '  ' ^  oh  the  union  of 
all  pairs  of n -constant  lines  which 
intersect  at  the  same  interior 
point  of  the  re-entrant  segment. 

(vi)  Using  a  local  coordinate  system,  the 

value  of  p',<’^  at  the  singular 
point  is  computed  explicitly. 

(vii)  Calculate  the  values  of  p^+1^ 

on  the  two n “constant  lines  which 
intersect  the  singjlar  point. 

Sc > era!  modifications  of  the  usual  AC!  scheme 
-i'i  necessary  due  to  the  recent-ant  boundary. 

A  periodic  trfdiajor.al  system,  was  solved  when 
ca-r.irg  out  the  ca'culation  in  (1).  The 
s.r'te"  in  (v)  contains  unknowns  on,  above, 
and  below  the  re-entrant  cut  in  the  physical 
region.  Therefore,  the  site  of  the  tridiagonal 
system  in  (v)  was  about  t-ice  as  large  as 
these  in  (iv)  and  (vii).  large  tridiagonal 
systems  of  this  type  would  be  cannon  in  the 
use  of  AC!  schemes  on  composite  grids.  The 
s,. stem,  in  general,  could  contain  unknowns 
from  several  different  subregions. 


In  this  particular  case,  a’l  the  unkn>nsof  a 
system  lie  on  n=constant  lines.  However,  the 
orientation  or  the  coordinate  line  does  change 
at  the  cut.  In  a  more  geneial  case,  some  of  the 
unknowns  could  follow  n*constant  lines  in'  one 
Subregion  and  C*constant  lines  in  another  sup- 
region.  df  course  these  same  concepts  would 
apply  to  any  of  the  other  fractional  step 
algorithms  such  as  approximate  factorization 
and  locally  one  dimensional  splitting. 

The  numerical ' solution  of  (10)  using  the 
iterative  algorithm  in  (12)  and  (13)  was 
compared  with  the  exact  solution,  and  the 
absolute  error  in  the  Surface  potential  is 
plotted  in  Figure  4.  It  is  interesting  to 
note  that  there  was  not  a  loss  In  accuracy 
at  the  outer  boundary.  Therefore,  the  finite- 
difference  solution  was  not  subjected  to  the 
same  type  of  translation  that  was  observed 
with  the  fini te -volume  solution.  The  similarities 
in  the  behavior  of  the  numerical  solution  at 
the  two  vertices  of  the  ellipse  Illustrate 
the  successful  incorporation  of  a  local  coordi¬ 
nate  system  at  the  singularity.  There  was 
only  a  slightly  larger  error  at  the  singular 
vertex. 


Several  interpolation  schemes  for  computing 
on  overlapping  grids  have  already  been  discussed. 
The  objective  new  is  to  compare  the  accuracy 
of  these  schemes.  The  following  example  also 
compares  the  solutions  for  overlapping  grids 
with  one  computed  on  a  composite  grid  with 
continuously  joined  grids.  The  function 


u(x,  y)  *  1  - 


log  (x2vy?) 
log  (200) 


satisfies  Laplace’s  equation  everywhere  except 
at  the  origin.  This  function  was  approximated 
by  solving  Laplace's  equation  on  the  region 
covered  by  the  overlapping  grid  in  Figure  5. 

The  exact  boundary  values  were  used.  Since 
the  Interior  boundary  component  was  taken  to 
be  the  unit  circle,  the  solution  assumed  the 
value  1  on  the  interior  boundary  and  decreased 
to  nearly  zero  at  the  outer  boundary.  The 
maximum  absolute  error  in  the  numerical 
solution  is  contained  in  Table  1.  Since  the 
linear  and  bilinear  interpolation  formulas 
can  be  applied  with  a  smaller  overlap  region, 
the  error  in  using  these  formulas  with  a 
one  cell  overlap  is  also  included.  This  problem 
was  also  solved  on  the  grid  in  Figure  6  which 
has  approximately  the  same  size  grid  cells.  • 

The  error  In  that  solution  is  also  Included 
in  Table  1.  All  the  computed  solutions  had 
about  the  same  order  of  accuracy.  The  second- 
order  interpolation  formulas  gave  better  results 
in  each  case,  although  the  difference  in  the 
error  was  minimal  when  the  same  overlap  was 
used. 
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A  more  interesting  comparison  arises  from 
tne  function 


u(x,  y)  •  sin  ==SL=== 

\x^  +  y£ 

which  is  the  solution  of  tne  Poisson  equation: 


where  +  is  the  appropriate  source  term. 

The  growth  of  the  higher-order  derivatives  results 
in  a  larger  truncation  error  in  this  example. 

The  maximum  absolute  error  in  the  numerical 
solution  for  some  of  the  overlap  cases  in  the 
previous  example  Is  listed  in  Table  2.  For 
this  function  the  quadratic  Taylor  polynomial 
gave  poorer  results  than  the  linear  Taylor 
polynomial.  Thus  we  have  a  case  where  the 
numerical  approximation  of  second-order 
derivatives  was  unreliable.  The  biquadratic 
interpolation  formula,  which  does  not  require 
numerical  differentiation,  was  clearly  superior 
to  bilinear  interpolation. 

All  of  these  results  were  computed  using  a 
point  SOS  iterative  method.  The  convergence 
rate  was  considerably  greater  with  an  overlap 
of  two  grid  cells  than  with  a  one  grid  cell 
Overlap.  These  results  are  thus  consistent 
with  previously  reported  observations  on  the 
relation  between  the  extent  of  o.erlap  and  the 
convergence  and  accuracy  of  the  numerical 
solution. 7  a  slightly  faster  rate  of  convergence 
was  also  noted  when  using  those  interpolation 
formulas  with  positive  coefficients,  namely, 
linear  and  bilinear. 


Conclusions 

Most  of  the  currently  available  algorithms 
for  the  numerical  solution  of  partial  differential 
equations  can  be  implemented  on  composite  grid 
systems.  Local  coordinate  systems  may  be  re¬ 
quired  at  special  points  where  the  usual 
differencing  techniques  cannot  be  applied. 

Finite  volume  formulations  are  easier  to  derive 
on  composite  grids,  and  In  principle,  may, 
be  derived  for  partial  differential  equations 
of  all  types.  When  using  ADI-type  algorithms, 
the  main  problem  Is  to  arrange  the  order  of 
the  computational  steps  so  that  all  the 
unknowns  at  each  level  are  updated  simultaneously. 
Except  when  using  ADI  schemes,  the  overlapping 
of  grids  is  an  alternative  to  the  more  common 
grid  construction  procedure  where  grid  lines 
continue  smoothly  from  or.e  subregion  Into  the 
next.  The  accuracy  of  numerical  solutions 
computed  on  overlapping  grids  Is  dependent 
on  both  the  interpolation  formula  and  the 
extent  of  the  overlap.  Several  Interpolation 
formulas  were  analyzed.  In  the  solution 
of  model  problems,  the  correct  choice  of  the 
interpolation  formula  could  reduce  the  error 
by  a  factor  of  two. 
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The  most  difficult  task  in  applying 
overlapping  grids  Is  to  partition  each  sub- 
region  grid  Into  three  sets;  those  points  where 
the  difference  equation  is  applied,  the  points 
whpre  the  Interpolation  formula  Is  used,  and  any 
points  which  are  not  used  in  the  numerical 
solution  of  the  prcbl.em.  The  code  which  was 
written  to  compute  the  examples  of  this  report 
attempted  to  automate  this  procedure  as  far  as 
possible.  In  the  process  several  different 
geometric  configurations  were  considered.  The 
solution  of  Laplace's  equation  with  values  0 
and  1  on  different  boundary  components  was 
computed  to  test  the  overlap  algorithm.  Three 
of  those  configurations  are  included  since 
they  are  representative  of  cases  where  over¬ 
lapping  grids  would  be  appropriate.  Figure 
7  typifies  a  body  In  a  channel  or  tunnel. 

Figure  8  is  representative  of  problems  involving 
flow  about  multiple  bodies.  Figure  9  Indicates 
a  case  where  overlapping  grids  can  be  used  to 
avoid  the  large  grid  line  curvature  that  would 
result  if  a  single  grid  were  constructed. 
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Interpolation  Formula  Uverlfo  *  J^or 
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Linear  on  triangles 
Bilinear 

linear  on  triangles 
Bi  1  inear 

Linear  Taylor  poly. 
Quadratic  on  tri. 
Biquadratic 
uadratic  Taylor  poly. 


ontinuously  joined  su 


0.00716 

0.00295 

0.00303 

0.00l58 

0.00273 

0.00200 

0.00156 

0.00211 


Table  1?  Error  in  the  approximation  of  u(x,y) 
l-log(xz+y*)/log{200)  by  solving  Laplace's 
equation. 


r  r. _  ,,  No.  Cells  Max.  Absolute 

Interpolation  Fomula  0y.  Error 


Bilinear 

Linear  Taylor  poly. 
Biquadratic 
uadratic  Taylor  pol 


0.05515 

0.02614 

0.02756 

0.03/58 


I ^b^sbs1.:: 


Table  2.  -Error  in  the  approximation  of  u{x,y)  ■ 
sin {2i//x£+y*)  by  solving  Poisson's  equation. 
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Figure  i.  Grid  for  computing  ideal  flow  about  an 
elliptic  cylinder;  complete  grid  and  closeup  near 
the  singular  point. 
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Figure  1.  Examples  of  special  points  and  their 
corresponding  local  coordinate  systems. 


Figure  3.  Error  in  the  surface  potential  computed 
using  the  finite-volume  method. 
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Numerical  Grid  Generation:  Foundations  and  Applications 


"Numerical  Implementation" 


- relevant  excerpt  from  Chapter  IV - 


by 


J.F.  Thompson,  Z.U.A.  Warsi,  C.W.  Mastin 


this  section  by 


H.V.  McConnaughey 


.  Discrete  Representation  of  Derivatives 

Approximate  values  of  the  spatial  derivatives  of  a 
function  which  appear  in  the  transformed  equations  may  be 
found  at  a  given  point  in  terms  of  the  function's  value  at 
that  point  and  at  neighboring  points.  As .noted  earlier, 
with  the  problem  in  the  transformed  space,  only  uniform 
square  grids  need  be  considered,  hence  the  standard  forms 
for  difference  representation  of  derivatives  may  be  used. 
For  example,  in  two  dimensions  the  first,  second,  and  mixed 
partials  with  respect  to  the  curvilinear  coordinates  £  and  q1 
are  ordinarily  represented  at  an  interior  point  (i,j)  by 
finite  differences  or  finite-volume  expressions  which  con¬ 
tain  function  values  at  no  more  than  the  nine  points  shown 
below. 

i-1  i  i+1 

o  o  «  j+1 


This  centered,  nine-point  "computational  molecule"  is 
usually  preferred  because  of  the  associated  difference 
representations  which  are  symmetry-  preserving  and  second- 
order  accurate.  Examples  of  finite-difference  approxima¬ 
tions  of  this  type  are: 


(Vij 

"  7(fi+1 ,j  ~  fi-1 ,j) 

(11a) 

(Vij 

"  2’(fi,j+1  ’  fi,J-1 5 

(lib) 

(f«)ij 

*  fi+1,j  '  2fij  +  fi-1 ,j 

(12a) 
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(12b) 


(fnn)iJ  “  fi,J+1  '  2fij  + 

“  if(fi+l.j+l  "  +  ri-l  f  J-1 5  (13) 


Other  second-order  approximations  of  the  mixed  partial 
which  use  the  nine-point  molecule  are: 

7(fi+1J+1  "  fi+1,J  ‘  fi,j+1  +  2fij  ‘  . 

‘  fi,j-1  '  fi-1,j  +  fi-1,j-1)  (111) 

and 

|(fi+1,j  '  fi+1,J-1  +  fi,j+1  ‘  2fij 
+  fi,j-1  "  fi-1,J+1  +  fi-1J)  (15) 


It  is  clear  that  at  boundary  points,  where  at'  most 
first  partials  must  be  represented,  the  computational  mole¬ 
cule  cannot  be  centered  relative  to  the  direction  of  the 
coordinate  5a  which  is  constant  on  the  boundary  (see  diagram 
below) . 


n  +  2 


n  +  iU 


m-i  m  m  +  i 
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There  a  one-sided  difference  must  be  used  to  approximate 

f  .  The  second-order  formula  appropriate  for  the  boundary 
ia 

point  indicated  above  is 


^f|Ct^mn  *  ^  ^m,n+2  +  ^m,n+1  “  3fm,n^ 


Any  standard  text'  on  the  subject  of  finite-difference  meth¬ 
ods  will  provide  formulas  of  alternate  order  and/or  based  on 
other  computational  molecules. 

A  finite-volume  approach  uses  function  values  at 
grid-cell  centers  and  approximates  derivatives  at  a  cell 
center  by  line  (surface  in  3D)  integrals  about  the  cell 
boundary  which  are  equivalent  to  averages  over  the  cell.  In 
particular,  the  identity 


^?avg 


ijIf  * 


D  do 


(16) 


is  used,  where  V  is  the  volume  of  D.  Thus,  if  a  function  is 
assumed  constant  along  a  grid-cell  face,  it  is  a  simple 
matter  to  evaluate  the  line  integral  in  (16)  when  D  is  a 
grid  cell  in  transformed  space.  In  terms  of  the  two- 
dimensional  grid: 
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this  approach  gives 


<fs>u 


-  f 


(17a) 


(fn)ij 


f  .  -  f 
2 


(1 7b) 


With  an  edge  value  approximated  as  the  average  of  the  center 
values  of  the  two  cells  sharing  that  face,  e.g. 


i  +  1 » j 


*  fu> 


(18) 


the  values  given  by  (17)  are  equivalent  to  ordinary  central 
differences  (cf.  Eq.  (11))  and  hence  are  second-order  ac¬ 
curate.  The  first  partials  of  f  may  also  be  assumed  con¬ 
stant  along  each  cell  edge  in  order  to  derive  from  (16)  the 
following  approximations  of  second  and  mixed  partials  at  a 
cell  center: 


(fVij 

-  <f-  .  -  fi)  . 

5  i+l,j  *  i-i,  j 

(19a) 

2  2 

(fWij 

-  (f-)  1  -  (fr)  1 

.  5  5  i.J-l 

2  2 

(19b) 

^n^ij 

-  (f_)  ,  -  (f_)  , 

n  n  i4.J 

2  2 

(19c) 

^nn^ij 

-  (f_)  ,  -  (f.)  , 

n  i » j+—  n  i,J-l 

2  2 

(19d) 

Now,  however,  the  averaging  scheme  in  (18)  cannot  be  used  to 
approximate  edge  values  of  the  derivatives  without  '  going 
outside  the  nine-point  computational  molecule  shown  above. 
Instead,  a  second-order  accurate  representation  can  be  ob¬ 
tained  on  the  nine-point  molecule  using  a  forward  (backward) 
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assignment  for  the  center  value  of  a  function  and  a  backward 
(forward)  assignment  for  the  first  partial  on  a  given  side. 
There  are  four  possible  schemes  of  this  type.  One  uses 


to  evaluate  Jf(E,n)  at  all  cell  centers  according  to  (17), 
and  then  uses 


gi-l,i  "  ’  gi  j-1  "  SU  (where  g  “  h  or  V  (21) 

to  evaluate  the  second  and  mixed  partials  given  in  (19). 
This  method  is  equivalent  to  a  finite-difference  scheme 
which  approximates  first  partials  by  backward  differences  of 
the  function,  and  then  approximates  second  and  mixed  par¬ 
tials  by  forward  differences  of  the  first  partials. 
Consequently,  the  second  derivatives  which  result  are  equal 
to  those  given  in  Eq.  (12),  while  the  resulting  representa¬ 
tions  of  the  two  mixed  partials  are  unequal  and  only  first- 
order  accurate.  If  the  two  mixed  partials  , are  averaged, 
however,  the  second-order  expression  (15)  is  recovered. 
This  is  also  true  of  the  reverse  scheme: 


(22) 


Vl,j  "  giJ’  gi,1+l  ’  SU  '  (g  “  fE  or  V  (23) 
2  2  ■ 

Expressions  (12)  and  (1*0  are  similarly  recovered  from  the 
other  two  possibilities  (Eq.  20a,  21a,  22b,  and  23b  or  Eq. 
20b,  21b,  22a,  and  23a).  The  symmetry-preserving  form  (13) 
can  be  recovered  by  averaging  the  averaged  mixed  partial  ob¬ 
tained  in  ,one  of  the  first  two  schemes  mentioned  and  that 
obtained  in  one  of  the  remaining  two. 
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The  manner  in  which  boundary  conditions  are  treated 
in  a  finite  volume  approach  depends  on  the  type  of  condi¬ 
tions  imposed.  When  Dirichlet  conditions  are  prescribed,  it 
is  advantageous  to  treat  the  boundary  as  the  center  line 
(plane  in  three-dimensions)  of  a  row  of  cells  straddling  the 
boundary.  The  centers  of  these  cells  then  fall  on  the  phys¬ 
ical  boundary  where  the  function  values  are  known.  When 
Neumann  or  mixed  conditions  are  given,  however,  'the  boundary 
is  best  treated  as  coincident  with  cell  faces. 

Suppose,  for  example,'  that  boundary  condition  (9)  is 
to  be1  imposed  at  the  cell  edge  indicated  below. 


i -Vi  i  +  >/2 


The  edge  value  of  cannot  be  approximated  by  the 

usual  averaging  scheme  (illustrated  by  Eq.  (18))  since  there 
is  no  cell  center  at  n-J~1 .  It'  can,  however,  be  found  in 
terms  of  neighboring  cell-centered  function  values  by  using 
boundary  condition  (9)  in  connection  with  the 
forward/backward  scheme  used  to  approximate  second  deriva¬ 
tives  at  the  cell  centers. 

Considering  the  '  scheme  represented  by  Eq.  (20)  and 
(21),  the  values  of  f  along  the  cell  edges  shown  above  are: 

*  f  i_i  4  f  i"f,i 

i-i,j  1  i ,  j+— 

2  2 


*  undetermined 


x 


i 


It  fellows  from  Eq.  (17)  that  the  first  partials  of  f  at  the 
cell  center  are 


(£Vij  *  fij  '  fi-1,J»  (Vij  ",  fij  '  xi 

Eq.  (21a, b)  then  give  f,  and. f^  along  the  cell  edges  enclos¬ 
ing  (i,j)  in  terms  of  fiHJ,  +  fH(j- 
x^  and  x^ .  In  particular, 


,2 


fij  ’ 


fi-i  ,j 


(fn)  i 

n  i» j  — 


fij  - 


Substitution  of  these  expressions  into  boundary  condition 
(9)  then  determines  the  edge  value  x^  as 


i.J 


1  -  (a  -  6u/g22)  1 


•  tr-4|rtg-”(r1J  -s22^]) 

/g— 


In  this  way,  f,  and  hence  f-  and  f^,  are  found  on  all 
boundary-cell  edges  in  terms  of  cell-centered  values  of  f. 

The  finite-difference  and  finite-volume  techniques 
described  thus  far  are  appropriate  for  representing  ail 
derivatives  with  respect  to  the  curvilinear  coordinates, 
even  those  appearing  in  the  metric  quantities.  In  fact,  as 
it  is  shown  later  in  this  chapter  and  in  chapter  V,  the  met¬ 
ric  quantities  should  be  represented  numerically  even  when 
analytical  expressions  are  available.  One  might  have,  for 
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example, 


^Vijk  “  i(xi  +  1  ,j,k  '  xi-1 , j.k*  (24) 
3.  Special  Points 

Many  of  the  expressions  given  in  the  previous  section 
break  down  at1 so-called  "special  points"  in  the  field  where 
special  attention  is  required  in  the  approximation  of  deriv¬ 
atives.  These  points  commonly  arise  when  geometrically 
complicated  physical  domains  are  involved.  As  indicated  in 
Chapter  II,  special  points  can  occur  on  the  domain  boundary 
and  on  interfaces  between  subregions  of  a  composite  curvi¬ 
linear  coordinate  system.  They  may  be  recognized  in  physi¬ 
cal  space  as  those  interior  points  having  a  nonstandard 
number  of  immediate’ neighbors  or,  equivalently,  those  points 
which  are  vertices,  or  the  center,  of  a  cell  with  either1  a 
nonstandard  number  of  faces  or  a  vertex  shared  by  a  nonstan¬ 
dard  number  .of  other  cells.  (In  two  dimensional  domains, 
ordinary  interior  points  have  eight  immediate  neighbors 
[refer  to  figure  on  p.  141];  standard  two-dimensional  in¬ 
terior  grid  cells  have  four  sides  and  share  each  vertex  with 
three  other  cells  [see  diagram  on  p.  143].)  Boundary  points 
are  not  special -unless  they  are  vertex-centered  and  have  a 
nonstandard  number  of  immediate,  neighbors  (other  than  five 
in  two  dimensions  see  diagram  on  p.  142  for  an  ordinary 
boundary  point)  and  then  are  special  only  when  their  asso¬ 
ciated  boundary  conditions  contain  spatial  derivatives.1 
.Some  examples  of  special’  cell-centered  points  and  special1 
vertex-centered  points  are  shown  below. 
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When  a  finite-difference  formulation  is  used,  the 
usual  approach,  as  described  in  Section  2,  can  be  followed 
at  a  special  point  P  if  the  transformed  equations  and 
difference  approximations  at  that  point  are  rephrased  in 
terms  of  suitable  local  coordinates.  The'  local  'system  is 
chosen  so  as  to  orient  and  label  only  the  surrounding  points 
to  be  used  m  the  needed  difference  expressipns.  Choices 
appropriate  to  various  special  points  are  listed  in  Tables 
1 ,  2 ,  and  3 • 

.  The  difficulties  encountered  ,at  special  points  in  a 
finite-volume  approach  are  clearly  seen  by  considering  the 
image  in  the  transformed  plane.  The  first  pair  of  diagrams 
below,  for  example,  shows  that  at  centers  cf  cells  having 
the  usual  number  of  faces  but  sharing, a  vertex  with  a  non¬ 
standard  number  of  cells,  such  difficulties  amount  to  mere 
bookkeeping  complications  when  only  first  partials  must  be 
approximated.  Equations  (17)  and  (18)  still  apply,  but  the 
indices  must  be  defined  to  correctly  relate  the  cell  centers 
on  the  two  sides  of  an  interface. ,  The  following  diagrams 
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Interior  body  and  is 
transformed  to  an  end 
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(f)  Special  point  5  Is  common  to 
three  subregions  and  is  an 
edge  point  for  each. 


Diagram  of 
special  point 

Characterization  of 
associated  special  cell 

Poin*  in  local 
computational  molecule 

vm 

Cell  center  is 
equivalent  to  special 
point  6  in  category 

IV. 

(a)  Points  1-9,  or 

(b)  points  1,  2,  4-6„ 

8,  9  and  use  the 
corresponding 
antisymmetric,  2nd- 
order  difference  for 
the  cross  derivative. 

DC  w 

x\>{  6  ! 
X>rTT 

Cell  center  is 
equivalent  to  special 
point  5  in  category 

V. 

(a)  Points  3-10  and  1 
or  2,  or 

(b)  3-7,  9,  10  and  use 
corresponding 
antisymmetric,  2nd- 
order  difference  for 
cross  derivative. 

x  V  / 

r/^4 — L 
\ 

. 

Cell  center  is 
equivalent  to  special 
point  5  in  category 

VI. 

‘ 

same  as  IX. 

XI  .  / 

sC/X'y*? 

Cell  center  is 
equivalent  to  special 
point  6  in  category 

VII. 

(i)  At  special  vertex 

3:  Use  points  1-8  and 

9  or  10. 

(ii)  Treat  special 
vertex  9  the  same  as 
special  vertex  4  in 
category  IX. 

Table  3.  Special  vertex-centered  interior  points 
associated  with  subregions  joined  between  grid  lines. 
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also  Illustrate  the  breakdown  at  all  special  cell-centered 
points  of  the  previously-described  finite-volume  schemes  for 
approximating  second  and  mixed  partial  derivatives.  This  is 
because  the  forward/backward  orientation  of  the  coordinate 
system  in  one  segment  cannot  be  consistently  followed  across 
the  interface  adjacent  to,  or  intersecting,  the  special 
points.  The  second  pair  of  diagrams  displays  the  additional 
complication  associated  with  grid  cells  having  a  nonstandard 
number  of  edges.  Such  a  cell  can  occur  on  an  interface  be¬ 
tween  segments  of  a  composite  grid  which  are  joined  between 
grid  lines.  When  the  segments  are  transformed  to  their  re¬ 
spective  images,  the  separate  pieces  of  the  special  grid 
cell  cannot  be  joined  without  distorting  them.  It  is  thus 
unclear  how  to  evaluate  the  volume  and  the  outward  normals 
of  that  transformed  cell  in  order  to  use  identity  (16)  in 
the  transformed  plane.  Consequently,  at  special  points  of 
this  type  and  at  all  special  points  where  second  derivatives 
must  be  approximated,  the  governing  equations  are  best 
represented  locally  in  the  physical  plane. where  such  ambi¬ 
guities  do  not  exist. 
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Treatment  in  physical  space  Involves  approximation  of 
the  original  equations  by  means  of  identity  (16).  Thus,  for 
a  two-dimensional  N-sided  cell  of  area  A  with  cartesian  cen¬ 
troid  P  -  (Pi»p2^»  vertices  V^-Cvj.v^)  i*1,2,...,N,  and 

edges  s*  joining  V*  and  V**1  (V^-V1)  along  which  a  func¬ 
tion  f  and  its  first  partial  derivatives  are  constant,  this 
approach  gives 

'l  ■  A'’  j,  ’  -  v|) 

f^x  "  ill  '  V2J 

N  i 

fP  -  A_1  7  f3  (yl  .  vi+1) 

yy  ry  tvi  vi  ; 

where  the  superscripts  on  f  and  its  derivatives  indicate  the 
point  or  face  of  evaluation.  As  in  the  previous  section,  an 
obvious  way  to  approximate  fs  is  to  average  the  center  val¬ 
ues  of  the  two  cells  sharing  edge  s1.  This  same  averaging 
scheme  cannot  be  repeated  to  approximate  f^  ,  and  fs  , 
however,  without  rejecting  the  recommended  strategy  of 
avoiding  use  of  values  at  points  which  are  not  immediate 
neighbors  of  the  point  at  which  a  quantity  is  being  evalu¬ 
ated.  Instead,  we  propose  the  averaging  technique: 


+  r 


,i*i 
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where  the  vertex  values  are  obtained  by  applying  identity 
(16)  to  auxiliary  cells  formed  by  joining  the  midpoints  of 
the  edges  of  each  cell  to  the  cell  center.  To  make  this 
more  precise,  let  V  be  a  vertex  common  to  Q  cells  and  label 
the  cell  faces  emanating  from  V  as  ki  with  midpoints 

Mi-(a1i,1m|)  1-1,2 . Q. 

Then  if  Pi-( p^- , p^ )  is  the  center  of  the  cell  having  edges  k1 
and  k^"*”1 ,  and  if 

f  -  fpi  along  M1?1  and  P 1Mi+1 

the  first  partial  derivatives  of  f  at  V  may  be  approximated 
by 

f»  -  A'1  l  -  mi) 

X  i-1  *  < 

fJ.A-'  f  fpI(»}-m}*’) 

J  i*1 

where  A  is  the  area  of  the  2Q-faced  auxiliary  cell 
.  indicated  in  the  following  diagram. 


This  technique  Is  applicable  to  all  grid  cell  centers; 
however,  It  Is  recommended  for  use  only  at  points  where 
methods  developed  In  section  2  break  down,  since  the  differ¬ 
ence  representations  associated  with  those  methods  are 

simpler. 
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